home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 4: GNU Archives / Linux Cubed Series 4 - GNU Archives.iso / gnu / glibc-1.09 / glibc-1 / glibc-1.09.1 / sysdeps / generic / log.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-07-31  |  14.1 KB  |  487 lines
  1. /*
  2.  * Copyright (c) 1992, 1993
  3.  *    The Regents of the University of California.  All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. #ifndef lint
  35. static char sccsid[] = "@(#)log.c    8.2 (Berkeley) 11/30/93";
  36. #endif /* not lint */
  37.  
  38. #include <math.h>
  39. #include <errno.h>
  40.  
  41. #include "mathimpl.h"
  42.  
  43. /* Table-driven natural logarithm.
  44.  *
  45.  * This code was derived, with minor modifications, from:
  46.  *    Peter Tang, "Table-Driven Implementation of the
  47.  *    Logarithm in IEEE Floating-Point arithmetic." ACM Trans.
  48.  *    Math Software, vol 16. no 4, pp 378-400, Dec 1990).
  49.  *
  50.  * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256,
  51.  * where F = j/128 for j an integer in [0, 128].
  52.  *
  53.  * log(2^m) = log2_hi*m + log2_tail*m
  54.  * since m is an integer, the dominant term is exact.
  55.  * m has at most 10 digits (for subnormal numbers),
  56.  * and log2_hi has 11 trailing zero bits.
  57.  *
  58.  * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h
  59.  * logF_hi[] + 512 is exact.
  60.  *
  61.  * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ...
  62.  * the leading term is calculated to extra precision in two
  63.  * parts, the larger of which adds exactly to the dominant
  64.  * m and F terms.
  65.  * There are two cases:
  66.  *    1. when m, j are non-zero (m | j), use absolute
  67.  *       precision for the leading term.
  68.  *    2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1).
  69.  *       In this case, use a relative precision of 24 bits.
  70.  * (This is done differently in the original paper)
  71.  *
  72.  * Special cases:
  73.  *    0    return signalling -Inf
  74.  *    neg    return signalling NaN
  75.  *    +Inf    return +Inf
  76. */
  77.  
  78. #if defined(vax) || defined(tahoe)
  79. #define _IEEE        0
  80. #define TRUNC(x)    x = (double) (float) (x)
  81. #else
  82. #define _IEEE        1
  83. #define endian        (((*(int *) &one)) ? 1 : 0)
  84. #define TRUNC(x)    *(((int *) &x) + endian) &= 0xf8000000
  85. #define infnan(x)    0.0
  86. #endif
  87.  
  88. #define N 128
  89.  
  90. /* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
  91.  * Used for generation of extend precision logarithms.
  92.  * The constant 35184372088832 is 2^45, so the divide is exact.
  93.  * It ensures correct reading of logF_head, even for inaccurate
  94.  * decimal-to-binary conversion routines.  (Everybody gets the
  95.  * right answer for integers less than 2^53.)
  96.  * Values for log(F) were generated using error < 10^-57 absolute
  97.  * with the bc -l package.
  98. */
  99. static double    A1 =       .08333333333333178827;
  100. static double    A2 =       .01250000000377174923;
  101. static double    A3 =     .002232139987919447809;
  102. static double    A4 =    .0004348877777076145742;
  103.  
  104. static double logF_head[N+1] = {
  105.     0.,
  106.     .007782140442060381246,
  107.     .015504186535963526694,
  108.     .023167059281547608406,
  109.     .030771658666765233647,
  110.     .038318864302141264488,
  111.     .045809536031242714670,
  112.     .053244514518837604555,
  113.     .060624621816486978786,
  114.     .067950661908525944454,
  115.     .075223421237524235039,
  116.     .082443669210988446138,
  117.     .089612158689760690322,
  118.     .096729626458454731618,
  119.     .103796793681567578460,
  120.     .110814366340264314203,
  121.     .117783035656430001836,
  122.     .124703478501032805070,
  123.     .131576357788617315236,
  124.     .138402322859292326029,
  125.     .145182009844575077295,
  126.     .151916042025732167530,
  127.     .158605030176659056451,
  128.     .165249572895390883786,
  129.     .171850256926518341060,
  130.     .178407657472689606947,
  131.     .184922338493834104156,
  132.     .191394852999565046047,
  133.     .197825743329758552135,
  134.     .204215541428766300668,
  135.     .210564769107350002741,
  136.     .216873938300523150246,
  137.     .223143551314024080056,
  138.     .229374101064877322642,
  139.     .235566071312860003672,
  140.     .241719936886966024758,
  141.     .247836163904594286577,
  142.     .253915209980732470285,
  143.     .259957524436686071567,
  144.     .265963548496984003577,
  145.     .271933715484010463114,
  146.     .277868451003087102435,
  147.     .283768173130738432519,
  148.     .289633292582948342896,
  149.     .295464212893421063199,
  150.     .301261330578199704177,
  151.     .307025035294827830512,
  152.     .312755710004239517729,
  153.     .318453731118097493890,
  154.     .324119468654316733591,
  155.     .329753286372579168528,
  156.     .335355541920762334484,
  157.     .340926586970454081892,
  158.     .346466767346100823488,
  159.     .351976423156884266063,
  160.     .357455888922231679316,
  161.     .362905493689140712376,
  162.     .368325561158599157352,
  163.     .373716409793814818840,
  164.     .379078352934811846353,
  165.     .384411698910298582632,
  166.     .389716751140440464951,
  167.     .394993808240542421117,
  168.     .400243164127459749579,
  169.     .405465108107819105498,
  170.     .410659924985338875558,
  171.     .415827895143593195825,
  172.     .420969294644237379543,
  173.     .426084395310681429691,
  174.     .431173464818130014464,
  175.     .436236766774527495726,
  176.     .441274560805140936281,
  177.     .446287102628048160113,
  178.     .451274644139630254358,
  179.     .456237433481874177232,
  180.     .461175715122408291790,
  181.     .466089729924533457960,
  182.     .470979715219073113985,
  183.     .475845904869856894947,
  184.     .480688529345570714212,
  185.     .485507815781602403149,
  186.     .490303988045525329653,
  187.     .495077266798034543171,
  188.     .499827869556611403822,
  189.     .504556010751912253908,
  190.     .509261901790523552335,
  191.     .513945751101346104405,
  192.     .518607764208354637958,
  193.     .523248143765158602036,
  194.     .527867089620485785417,
  195.     .532464798869114019908,
  196.     .537041465897345915436,
  197.     .541597282432121573947,
  198.     .546132437597407260909,
  199.     .550647117952394182793,
  200.     .555141507540611200965,
  201.     .559615787935399566777,
  202.     .564070138285387656651,
  203.     .568504735352689749561,
  204.     .572919753562018740922,
  205.     .577315365035246941260,
  206.     .581691739635061821900,
  207.     .586049045003164792433,
  208.     .590387446602107957005,
  209.     .594707107746216934174,
  210.     .599008189645246602594,
  211.     .603290851438941899687,
  212.     .607555250224322662688,
  213.     .611801541106615331955,
  214.     .616029877215623855590,
  215.     .620240409751204424537,
  216.     .624433288012369303032,
  217.     .628608659422752680256,
  218.     .632766669570628437213,
  219.     .636907462236194987781,
  220.     .641031179420679109171,
  221.     .645137961373620782978,
  222.     .649227946625615004450,
  223.     .653301272011958644725,
  224.     .657358072709030238911,
  225.     .661398482245203922502,
  226.     .665422632544505177065,
  227.     .669430653942981734871,
  228.     .673422675212350441142,
  229.     .677398823590920073911,
  230.     .681359224807238206267,
  231.     .685304003098281100392,
  232.     .689233281238557538017,
  233.     .693147180560117703862
  234. };
  235.  
  236. static double logF_tail[N+1] = {
  237.     0.,
  238.     -.00000000000000543229938420049,
  239.      .00000000000000172745674997061,
  240.     -.00000000000001323017818229233,
  241.     -.00000000000001154527628289872,
  242.     -.00000000000000466529469958300,
  243.      .00000000000005148849572685810,
  244.     -.00000000000002532168943117445,
  245.     -.00000000000005213620639136504,
  246.     -.00000000000001819506003016881,
  247.      .00000000000006329065958724544,
  248.      .00000000000008614512936087814,
  249.     -.00000000000007355770219435028,
  250.      .00000000000009638067658552277,
  251.      .00000000000007598636597194141,
  252.      .00000000000002579999128306990,
  253.     -.00000000000004654729747598444,
  254.     -.00000000000007556920687451336,
  255.      .00000000000010195735223708472,
  256.     -.00000000000017319034406422306,
  257.     -.00000000000007718001336828098,
  258.      .00000000000010980754099855238,
  259.     -.00000000000002047235780046195,
  260.     -.00000000000008372091099235912,
  261.      .00000000000014088127937111135,
  262.      .00000000000012869017157588257,
  263.      .00000000000017788850778198106,
  264.      .00000000000006440856150696891,
  265.      .00000000000016132822667240822,
  266.     -.00000000000007540916511956188,
  267.     -.00000000000000036507188831790,
  268.      .00000000000009120937249914984,
  269.      .00000000000018567570959796010,
  270.     -.00000000000003149265065191483,
  271.     -.00000000000009309459495196889,
  272.      .00000000000017914338601329117,
  273.     -.00000000000001302979717330866,
  274.      .00000000000023097385217586939,
  275.      .00000000000023999540484211737,
  276.      .00000000000015393776174455408,
  277.     -.00000000000036870428315837678,
  278.      .00000000000036920375082080089,
  279.     -.00000000000009383417223663699,
  280.      .00000000000009433398189512690,
  281.      .00000000000041481318704258568,
  282.     -.00000000000003792316480209314,
  283.      .00000000000008403156304792424,
  284.     -.00000000000034262934348285429,
  285.      .00000000000043712191957429145,
  286.     -.00000000000010475750058776541,
  287.     -.00000000000011118671389559323,
  288.      .00000000000037549577257259853,
  289.      .00000000000013912841212197565,
  290.      .00000000000010775743037572640,
  291.      .00000000000029391859187648000,
  292.     -.00000000000042790509060060774,
  293.      .00000000000022774076114039555,
  294.      .00000000000010849569622967912,
  295.     -.00000000000023073801945705758,
  296.      .00000000000015761203773969435,
  297.      .00000000000003345710269544082,
  298.     -.00000000000041525158063436123,
  299.      .00000000000032655698896907146,
  300.     -.00000000000044704265010452446,
  301.      .00000000000034527647952039772,
  302.     -.00000000000007048962392109746,
  303.      .00000000000011776978751369214,
  304.     -.00000000000010774341461609578,
  305.      .00000000000021863343293215910,
  306.      .00000000000024132639491333131,
  307.      .00000000000039057462209830700,
  308.     -.00000000000026570679203560751,
  309.      .00000000000037135141919592021,
  310.     -.00000000000017166921336082431,
  311.     -.00000000000028658285157914353,
  312.     -.00000000000023812542263446809,
  313.      .00000000000006576659768580062,
  314.     -.00000000000028210143846181267,
  315.      .00000000000010701931762114254,
  316.      .00000000000018119346366441110,
  317.      .00000000000009840465278232627,
  318.     -.00000000000033149150282752542,
  319.     -.00000000000018302857356041668,
  320.     -.00000000000016207400156744949,
  321.      .00000000000048303314949553201,
  322.     -.00000000000071560553172382115,
  323.      .00000000000088821239518571855,
  324.     -.00000000000030900580513238244,
  325.     -.00000000000061076551972851496,
  326.      .00000000000035659969663347830,
  327.      .00000000000035782396591276383,
  328.     -.00000000000046226087001544578,
  329.      .00000000000062279762917225156,
  330.      .00000000000072838947272065741,
  331.      .00000000000026809646615211673,
  332.     -.00000000000010960825046059278,
  333.      .00000000000002311949383800537,
  334.     -.00000000000058469058005299247,
  335.     -.00000000000002103748251144494,
  336.     -.00000000000023323182945587408,
  337.     -.00000000000042333694288141916,
  338.     -.00000000000043933937969737844,
  339.      .00000000000041341647073835565,
  340.      .00000000000006841763641591466,
  341.      .00000000000047585534004430641,
  342.      .00000000000083679678674757695,
  343.     -.00000000000085763734646658640,
  344.      .00000000000021913281229340092,
  345.     -.00000000000062242842536431148,
  346.     -.00000000000010983594325438430,
  347.      .00000000000065310431377633651,
  348.     -.00000000000047580199021710769,
  349.     -.00000000000037854251265457040,
  350.      .00000000000040939233218678664,
  351.      .00000000000087424383914858291,
  352.      .00000000000025218188456842882,
  353.     -.00000000000003608131360422557,
  354.     -.00000000000050518555924280902,
  355.      .00000000000078699403323355317,
  356.     -.00000000000067020876961949060,
  357.      .00000000000016108575753932458,
  358.      .00000000000058527188436251509,
  359.     -.00000000000035246757297904791,
  360.     -.00000000000018372084495629058,
  361.      .00000000000088606689813494916,
  362.      .00000000000066486268071468700,
  363.      .00000000000063831615170646519,
  364.      .00000000000025144230728376072,
  365.     -.00000000000017239444525614834
  366. };
  367.  
  368. double
  369. #ifdef _ANSI_SOURCE
  370. log(double x)
  371. #else
  372. log(x) double x;
  373. #endif
  374. {
  375.     int m, j;
  376.     double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0;
  377.     volatile double u1;
  378.  
  379.     /* Catch special cases */
  380.     if (x <= 0)
  381.         if (_IEEE && x == zero)    /* log(0) = -Inf */
  382.             return (-one/zero);
  383.         else if (_IEEE)        /* log(neg) = NaN */
  384.             return (zero/zero);
  385.         else if (x == zero)    /* NOT REACHED IF _IEEE */
  386.             return (infnan(-ERANGE));
  387.         else
  388.             return (infnan(EDOM));
  389.     else if (!finite(x))
  390.         if (_IEEE)        /* x = NaN, Inf */
  391.             return (x+x);
  392.         else
  393.             return (infnan(ERANGE));
  394.     
  395.     /* Argument reduction: 1 <= g < 2; x/2^m = g;    */
  396.     /* y = F*(1 + f/F) for |f| <= 2^-8        */
  397.  
  398.     m = logb(x);
  399.     g = ldexp(x, -m);
  400.     if (_IEEE && m == -1022) {
  401.         j = logb(g), m += j;
  402.         g = ldexp(g, -j);
  403.     }
  404.     j = N*(g-1) + .5;
  405.     F = (1.0/N) * j + 1;    /* F*128 is an integer in [128, 512] */
  406.     f = g - F;
  407.  
  408.     /* Approximate expansion for log(1+f/F) ~= u + q */
  409.     g = 1/(2*F+f);
  410.     u = 2*f*g;
  411.     v = u*u;
  412.     q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
  413.  
  414.     /* case 1: u1 = u rounded to 2^-43 absolute.  Since u < 2^-8,
  415.      *            u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits.
  416.      *         It also adds exactly to |m*log2_hi + log_F_head[j] | < 750
  417.     */
  418.     if (m | j)
  419.         u1 = u + 513, u1 -= 513;
  420.  
  421.     /* case 2:    |1-x| < 1/256. The m- and j- dependent terms are zero;
  422.      *         u1 = u to 24 bits.
  423.     */
  424.     else
  425.         u1 = u, TRUNC(u1);
  426.     u2 = (2.0*(f - F*u1) - u1*f) * g;
  427.             /* u1 + u2 = 2f/(2F+f) to extra precision.    */
  428.  
  429.     /* log(x) = log(2^m*F*(1+f/F)) =                */
  430.     /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q);    */
  431.     /* (exact) + (tiny)                        */
  432.  
  433.     u1 += m*logF_head[N] + logF_head[j];        /* exact */
  434.     u2 = (u2 + logF_tail[j]) + q;            /* tiny */
  435.     u2 += logF_tail[N]*m;
  436.     return (u1 + u2);
  437. }
  438.  
  439. /*
  440.  * Extra precision variant, returning struct {double a, b;};
  441.  * log(x) = a+b to 63 bits, with a is rounded to 26 bits.
  442.  */
  443. struct Double
  444. #ifdef _ANSI_SOURCE
  445. __log__D(double x)
  446. #else
  447. __log__D(x) double x;
  448. #endif
  449. {
  450.     int m, j;
  451.     double F, f, g, q, u, v, u2, one = 1.0;
  452.     volatile double u1;
  453.     struct Double r;
  454.  
  455.     /* Argument reduction: 1 <= g < 2; x/2^m = g;    */
  456.     /* y = F*(1 + f/F) for |f| <= 2^-8        */
  457.  
  458.     m = logb(x);
  459.     g = ldexp(x, -m);
  460.     if (_IEEE && m == -1022) {
  461.         j = logb(g), m += j;
  462.         g = ldexp(g, -j);
  463.     }
  464.     j = N*(g-1) + .5;
  465.     F = (1.0/N) * j + 1;
  466.     f = g - F;
  467.  
  468.     g = 1/(2*F+f);
  469.     u = 2*f*g;
  470.     v = u*u;
  471.     q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
  472.     if (m | j)
  473.         u1 = u + 513, u1 -= 513;
  474.     else
  475.         u1 = u, TRUNC(u1);
  476.     u2 = (2.0*(f - F*u1) - u1*f) * g;
  477.  
  478.     u1 += m*logF_head[N] + logF_head[j];
  479.  
  480.     u2 +=  logF_tail[j]; u2 += q;
  481.     u2 += logF_tail[N]*m;
  482.     r.a = u1 + u2;            /* Only difference is here */
  483.     TRUNC(r.a);
  484.     r.b = (u1 - r.a) + u2;
  485.     return (r);
  486. }
  487.